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INTRODUCTION 


The  present  document  is  a  final  report  for  work  on  supersonic 
nozzle  analysis  for  both  perfect  and  real  gases  and  for  viscous  and 
inviscid  flows.  The  primary  focus  of  this  research  was  to  develop 
an  accurate,  efficient  procedure  for  computing  nozzle  flows  with 
careful  attention  to  global  mass  conservation.  In  general,  the 
goal  was  to  demonstrate  computations  with  mass  flux  errors  of  less 
than  1%  on  moderately  resolved  grids.  Calculation  procedures  for 
both  viscous  and  inviscid  flows  and  for  parabo] ized  Navier-Stokes 
( PNS )  algorithms  have  consistently  demonstrated  this  goal.  Primary 
emphasis  has  been  on  perfect  gas  calculations,  but  a  real  gas 
analysis  has  been  developed  and  demonstrated  to  a  lesser  degree. 

Also  included  in  the  effort  was  a  procedure  for  coupling  the  nozzle 
wall  cooling  flow  with  the  nozzle  flowfield  so  that  both  nozzle 
wall  temperature  and  heat  flux  were  simultaneously  determined  by 
the  calculation.  Finally,  a  procedure  for  determining  the  effect 
of  small  variations  in  back  pressure  on  the  nozzle  flowfield  has 
been  developed. 

The  present  final  report  includes  as  Appendix  C  a  Ph.D. 
thesis*  which  describes  the  detail  of  the  perfect  gas  work.  In 
particular,  this  Appendix  includes  the  coupled  nozzle  wall  cooling 
and  nozzle  back  pressure  calculations  as  an  appendix.  Results  of 
this  work  have  also  appeared  as  parts  of  three  papers2/ 3/ 4  which 
are  available  in  the  literature.  The  main  body  cf  the  report 
contains  nozzle  back  pressure  calculations,  including  flows  with 
weakly  separated  boundary  layers  and  associated  recirculation 
regions.  This  Appendix  also  describes  methods  for  extending  the 
present  results  to  the  fully  three-dimensional  problem. 

The  main  body  of  the  report  is  focussed  on  the  real  gas 
formulation  and  on  representative  results  for  real  gas  calculations. 
In  addition,  some  recent  perfect  gas  calculations  which  do  not 
appear  in  Appendix  C  are  included. 

The  high  temperatures  associated  with  rocket  combustion 
processes  lead  to  vibrational  excitation  and  dissociation  of  the 
flowing  gases,  and  as  the  gas  expands  through  the  nozzle,  these 
high  temperature  effects  have  a  significant  influence  on  the 
physics  of  the  flow  and,  therefore,  cannot  be  computed  by  perfect 
gas  assumptions. 

The  most  complete  model  for  describing  these  phenomena  is  to 
add  the  appropriate  species  diffusion  equations  and  take  into 
account  the  complete  chemical  kinetic  effects.  The  present 
analysis,  however,  addresses  the  simpler  problem  of  a  real  gas 
undergoing  shifting  equilibrium  chemistry.  This  implies  that  the 
real  gas  chemistry  is  infinitely  rapid.  During  the  expansion,  the 
gases  in  the  nozzle  generally  reach  a  point  at  which  their 
temperature  and  pressure  are  low  enough  that  this  infinite  rate 
assumption  fails  and  non-equilibrium  effects  begin  to  appear.  The 
simplest  expedient  for  handling  this  phenomena  is  to  switch  from 
equilibrium  to  frozen  flow  assumptions  at  a  given  point.  Such 
issues  are,  however,  not  included  in  the  present  analysis  which 
discusses  only  the  equilibrium,  real  gas  problem. 

In  the  next  section,  the  governing  systems  of  equations  are 
described,  and  the  real  gas  model  is  defined.  The  real  gas  effects 
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assume  that  the  gas  properties,  including  the  molecular  weight,  Mw, 
the  temperature,  T,  the  internal  energy,  e,  the  viscosity,  p,  and 
the  thermal  conductivity,  K,  are  functions  of  the  thermodynamic 
variables.  For  completeness,  both  the  Euler  equations  and  the 
thin-layer  Navier-Stokes  (TLNS)  equations  are  considered. 

In  the  following  section,  the  implicit  numerical  schemes  for 
solving  the  Euler  and  TLNS  equations  are  described.  The  Euler 
solutions  are  limited  to  fully  supersonic  flows,  while  the  TLNS 
solutions  contain  a  thin  subsonic  layer  along  the  wall. 
Corresponding  transonic  codes  (as  descrioed  in  Appendix  C)  are  used 
to  generate  the  supersonic  (or  predominantly  supersonic)  start 
lines  for  these  supersonic-oriented  algorithms. 

An  important  difference  between  the  perfect  gas  analysis  of 
Appendix  C  and  the  real  gas  equations  of  the  main  body  of  the 
report  are  that  the  flux  vectors  of  the  former  are  homogeneous 
while  those  for  the  real  gas  equations  are  non-homogeneous .  This 
requires  that  a  different  splitting  of  the  Jacobians  associated 
with  convective  terms  be  used  for  the  difference  equations. 

Differential  equations  to  compute  the  nozzle  rate  of  mass  flow 
and  thrust  are  also  evaluated  in  this  section.  The  next  section 
includes  the  user's  manual  which  explains  how  to  run  the  supersonic 
codes.  The  following  section  presents  representative  real  gas 
results.  For  simplicity,  these  calculations  are  based  on  the 
properties  of  equilibrium  air  which  were  more  readily  available 
than  combustion  gas  properties.  The  extension  to  combustion  gas 
properties  is  straightforward  once  the  thermodynamic  properties  are 
available  for  the  desired  fuel,  oxidizer  and  equivalence  ratio. 
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THE  GOVERNING  EQUATIONS 


The  Euler  and  TINS  Equations 


The  two-dimensional  axisymmetric  form  of  the  unsteady 
Navier-Stokes  equations  written  in  generalized  body- fitted 
transformed  coordinates  are  given  by: 

.  9e  ,  3f  ,  3r  ,  3s 

at  +  ^  +  m  =  H  +  ^  + 
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Here,  (x,y)  are  the  Cartesian  coordinates  parallel  and 
perpendicular  to  the  axis  of  symmetry,  respectively;  (£,r|)  are  the 
body-fitted  transformed  coordinates, 


S  =  S(x,y) 
T]  =  n(*,y) 


(2) 


where  £  is  the  streamwise  coordinate  which  coalesces  with  the 
nozzle  axis  and  the  nozzle  wall;  J  is  the  transformation  Jacobian, 
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(3) 


3  *  xZy1  '  Vt 

and  U  and  V  are  the  contravariant  velocities, 

U  =  Zxu  +  £yV  and  V  =  iixu  +  HyV  (4) 

The  viscous  and  conductive  terms  are, 

’xx  “  P[2(Ex  at  +  lx  ly>  '  f  v'^] 

Tyx  “  Txy  =  31  +  ’’x  37>  +  <^y  at  +  "y  3n>  1 

ryy  -  Pt2^y  31  +  9y  3^>  '  3  ,5) 

3rp  Orp 
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2rn  3ip 
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v'y’  =  <E;x  at  +  lx  at*  +  y  +  ay  at  +  "y  3n> 

The  vector,  Q,  contains  the  dependent  variables.  Here,  the 
standard  notation  has  been  used  for  density,  p,  velocity 
components,  (u,v) ,  pressure,  p,  temperature,  T,  and  total  in  .arnal 
energy,  eQ.  The  vectors  E  and  F  are  the  convective  fluxes,  R  and  S 
are  the  viscous  fluxes,  and  H  is  the  source  term. 

The  equation  of  state  is, 

p  -  p  r  T  <6> 

w 

This  is  the  ideal  gas  equation  where  R  is  the  universal  gas 
constant. 

The  internal  energy  is  defined  in  terms  of  the  total  energy, 
the  density,  and  the  velocity  components, 
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K  =  K  (p,  e) 


conductivity 


(11) 


The  gas  properties  can  be  represented  in  functional  or  in  tabulated 
data  form.  For  a  perfect  gas,  Mw  is  constant  and  the  temperature 
can  be  related  to  the  other  variables  by, 


T  =  h  (f° '  2  (u2  +  v2)  (12) 

where  Cp  is  the  specific  heat  at  constant  pressure  and  Y  is  the 
specific  heat  ratio. 

The  Euler  equations  are  derived  from  the  full  Navier-Stokes 
equations  by  eliminating  all  the  viscous  and  heat  transfer  terms, 
i.e.,  p  =  k  =  0.  The  inviscid  equations  then  become  identical  to 
Eqn.  1  with  R  and  S  set  to  zero. 

The  advantage  of  this  reduced  set  of  equations  is  that  it  can 
be  solved  numerically  using  much  less  computer  time  and  memory  than 
is  required  for  the  complete  Navier-Stokes  equations. 

The  thin-layer  Navier-Stokes  (TLNS)  equations  take  into 
account  only  those  viscous  and  conductive  terms  containing 
derivatives  in  the  direction  perpendicular  to  the  body  surface. 
Terms  containing  derivatives  parallel  to  the  surface  are  dropped. 
Specifically,  the  thin-layer  approach  assumes  that  9/9£  =  0.  The 
TLNS  equations  are  thus  reduced  to. 


cto  9e  9f  _  ds' 

3t  +  3f  +  3^-H  +  g^T 


3s  ‘ 
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where  the  vectors  H'  and  S'  are  given  by, 
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where, 
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The  viscous  term,  S',  can  be  further  rearranged  as, 

=  3  ctoi  3^2 

3iT  Si]  (R1  Sn  2  sr} 


in  which  R^  and  R2  are  4x4  matrices  given  by, 


The  vectors,  Qi  and  Q2 ,  are  defined  by, 

Ql  =  (P/  u,  v,  eQ)T,  Q2  =  (T,  u2,  v2,  uv)T  (17) 

In  this  form,  the  viscous  dissipation  in  the  energy  equation  is 
separated  from  the  remaining  viscous  terms  and  the  matrices  Rj  and 
R2  contain  only  the  viscosity  and  thermal  conductivity  and  the 
metrics  of  the  transformation. 

Important  Relations  For  The  Real  Gas  Equations 

The  boundary  layers  in  rocket  nozzles  can  be  either  laminar  or 
turbulent,  and  frequently  it  is  difficult  to  tell  what  the  nature 
of  the  flow  will  be.  The  favorable  pressure  gradient,  the  wall 
cooling  and  the  supersonic  character  of  the  flow  all  tend  to 
stabilize  the  boundary  layer  and  favor  laminar  flow,  but  for  high 
Reynolds  numbers  (such  as  occur  in  large  thrust  systems) ,  and  high 
expansion  ratio  nozzles  (which  are  typically  quite  long) ,  the 
boundary  layers  are  generally  turbulent.  Modeling  turbulent  flows 
is  never  an  easy  assignment,  but  for  these  supersonic,  strongly 
accelerating,  highly  cooled  boundary  layers  with  chemical  reaction 
and  heat  release,  it  is  safe  to  speculate  that  existing  turbulence 
models  which  have  been  developed  and  calibrated  for  incompressible 
flows  with  weak  pressure  gradients  will  be  less  realistic  than  in 
their  normal,  more  benign,  environments. 

Consequently,  the  use  of  complex  turbulence  models  which  are 
interested  in  the  present  regime  is  hardly  justifiable. 

Accordingly,  we  have  used  the  simple  Baldwin-Lomax  two-layer  mixing 
length  model,  along  with  a  constant  turbulent  Prandtl  number 
assumption.  This  approach  gives  flowfield  profiles  which  are 
representative  of  turbulence  conditions  and  verifies  that  the  codes 
continue  to  work  in  the  turbulent  flow  environment.  Incorporation 
of  more  advances  turbulence  models  can  readily  be  incorporated  if 
and  when  they  become  available. 
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NUMERICAL  SOLUTION  PROCEDURES 


Although  the  fluid  dynamic  equations  are  formulated  in  the 
complete  unsteady  sense,  only  steady  state  solutions  are  considered. 
These  solutions  are  obtained  by  applying  time-dependent  schemes  to 
the  Euler  and  TLNS  equations.  Implicit  time-marching  methods  are 
used  throughout. 

The  spatial  discretization  that  is  chosen  for  the  axisymmetric 
problem  uses  central  differencing  in  the  cross-stream  or  radial  (i]) 
direction,  and  second  order  upwind  differencing  in  the  streamwise 
or  axial  (£)  direction.  For  the  inviscid  supersonic  case,  this 
upwind-central  differencing  combination  nicely  mimics  the  physics 
of  the  flow.  For  the  parabolized  case,  it  also  mimics  the  type  of 
differencing  that  is  typically  chosen  for  PNS  algorithms.  The 
extension  of  this  hybrid  differencing  scheme  to  the  TLNS  equations 
is  also  straightforward  if  we  interpret  the  upwind  differencing  in 
a  flux  split  sense. 

With  this  choice  of  differencing,  both  the  supersonic  inviscid 
case  and  the  PNS  formulation  can  be  solved  by  a  single  sweep 
through  the  flowfield.  The  TLNS  formulation  then  becomes  a  fully 
iterative  (forward-backward  sweeps)  solution  procedure  for  which 
the  PNS  solution  is  just  the  forward  sweep. 

Details  of  this  discretization  and  solution  method  are  given 
in  the  present  section  fdor  the  real  gas  formulation.  Similar 
developments  for  the  perfect  gas  case  are  given  in  Appendix  D. 
Because  the  flux  vectors  in  the  perfect  gas  formulation  are 
homogeneous  functions  of  the  dependent  variable,  Q,  flux  vector 
splitting,  rather  than  flux  difference  splitting,  is  used  in 
Appendix  D.  The  flux  difference  method  described  here  is  necessary 
for  real  gases  where  the  flux  vectors  are  non-homogeneous,  but 
should  perform  as  well  as  flux  vector  splitting  in  the  perfect  gas 
case.  Consequently,  our  recommendation  is  to  use  flux-difference 
splitting  throughout  rather  than  only  in  the  real  gas  case. 

We  begin  by  considering  the  thin-layer  Navier-Stokes  equations 
and  then  drop  viscous  effects  to  obtain  an  Euler  equation 
procedure. 


Discretization  Of  The  Thin-Layer  Navier-Stokes  Equations 


We  begin  by  discretizing  the  TLNS  equations  (Eqn.  12)  in  time 
using  Euler  implicit  differencing  in  time, 
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n+1 


=  0 


(18) 


where  n  and  n+1  refer  to  the  previous  and  current  time  levels, 
respectively.  We  then  express  Eqn.  18  in  delta  form, 

AQ  +  At(^AE  +  j^AF  -  ^  AS'  -  AH  =  -At(||  +  “  H')n  <19> 

where  AQ  =  Qn+1  -  Qn,  etc. 


The  spatial  discretization  of  the  right-hand  side  of  Eqn. 
19  is  accomplished  as  follows, 


9e  =  El+1/2^  -  Ei-1/?j 


§1  =  ^ij  +  1/2  ~  Fjj-1/?. 

oil  Ati 


and, 


-  <ri  . . 

ij+1/2  x  ai>  i j -1/2 


]i 


Ati 


(20) 


where  the  tilde  denotes  the  numerical  flux  of  the  indicated  vectors 
These  quantities  are  identified  below. 

For  the  streamwise  vector,  E,  the  discretization  is  defined 
by. 


Ei+l/2 j  2  tEij+  Ei+lj"  AEi+l/2j+  AEi+l/2j+  AEi-l/2j 

whereby  a  Taylors  series  in  £,  we  define. 


iEi+3/2j|21) 


AEi+l/2j  A  <Q)i+l/2j  ^i+lj  " 


where  A+  and  A“  are  the  eigenvalue  split  forms  of  A.  When  Eqn.  21 
is  substituted  into  the  discretized  form  of  3e/8£  defined  in  Eqn. 
20,  this  gives  a  flux  split  second  order  upwind  differenced  scheme. 

The  splitting  of  matrix  A  is  determined  on  the  basis  of  the 
eigenvalues  of  A.  Thus,  we  define, 


A  =  A+  +  A"  (22a) 

where  the  eigenvalues  of  A+  are  positive  and  those  of  A-  are 
negative.  We  define  this  matrix  by  diagonalizing  A  by  the 
similarity  transformation, 


Aa  =  T"1P~1APT 


where , 


Aa  =  diag  (U,  U,  U+c,  U-C) 


(23) 


Here,  U  is  the  contravariant  velocity,  U  £yv,  and  C  is  the 


2 


+  Sv2) 


where  c  is  the 


transformed  speed  of  sound,  C  =  c( Ex¬ 
physical  speed  of  sound. 

The  split  eigenvalues  are  defined  by  A^  =  ( A i  ±  Ai)/2  where 
Ai  are  the  individual  eigenvalues  in  Eqn.  23.  As  can  be  seen  by 
inspection,  this  choice  leads  to  the  definition  of  two  diagonal 
eigenvalues,  AA+  and  AA”  where  all  eigenvalues  of  AA+  are  either 
positive  or  zero  and  those  of  AA“  are  negative  or  zero.  From  these 
split  diagonal  matrices,  we  can  define  the  splitting  of  A  as. 


A+  =  TPAA+  P-1T-1 

A"  =  TPAa"  P-1T_1  (24) 


and  by  inspection  we  see  A+  +  A~  =  A.  The  matrix  PT  is  given  in 
Appendix  A. 
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Thus  far,  we  have  demonstrated  the  method  of  splitting  A  and 
its  upwind  differencing.  We  now  return  to  Eqn.  20  to  define  the 
spatial  differencing  in  rj .  Because  we  use  central  differencing  in 
H ,  we  define. 


Fij+l/2  =  l/2(Fij+i  +  Fij) 


which  leads  to  the  standard  three  point  differencing, 
differencing  of  the  two  viscous  term  is  given  as, 

h  Ri  s?1  = 


The 


i, .  <Sij^  Qi*>  -  Ri, 


-ij+1/2 


ij-1/2 


(25) 


(26) 


The  differencing  of  the  first  order  derivations  in  H'  are  done 
in  a  manner  identical  to  that  for  9f/9i).  The  derivatives  on  the 
left-hand  side  of  Eqn.  19  are  accomplished  with  the  aid  of  a 
Taylor's  series  in  time.  For  the  vector,  E,  we  have, 


En+1 


En  + 


9e 

3q  aq 


so  that, 

AE  =  En+1  -  En  =  ||  AQ 


(27) 

(28) 


Similar  results  can  be  obtained  for  the  vectors  F,  ,  Q2,  and  H'. 
The  notation  for  the  Jacobians  of  E  and  F  has  already  been 
introduced  (as  A  and  B,  respectively) .  We  also  define  the 
Jacobians, 


(29) 


These  matrices  are  given  in  Appendix  A.  We  also  use  the  flux  split 
matrices  A+  and  A”  for  the  streamwise  derivatives.  Incorporating 
these  Taylor's  series  expansions  on  the  left-hand  side  and  the 
previously  defined  differences  for  the  right-hand  side,  the  fully 
discretized  equation  becomes, 

[I  =  At(^(A++A~)  +  ^B  -  Rx  ^  Bvl  -  ^  R2  ^  Bv2 

-  D]  =  -At[||  +  ^  Rx  ^  R2  ~  H']  (30) 

For  simplicity  of  notation,  we  express, 

D  =  D'  +  D"  (31) 


where  D'  represents  the  Jacobian  of  the  algebraic  terms  in  H'  (Eqn. 
13)  and  D''  represents  the  Jacobian  of  the  first  derivative  terms 
in  H'.  These  matrices  are  defined  in  Appendix  A.  We  then  also 
express  the  T)  derivatives  on  the  left-hand  side  of  Eqn.  30  by  B<r 
(for  total  derivative)  so  that, 
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(32) 


h  btaq  =  Lh(B "  Ri  h  Bvi '  r2  h  bv2) 


D' '  ]AQ 


and  again  use  the  S'  notation  for  the  viscous  terms  on  the 
right-hand  side  of  Egn.  30  (see  Egn.  15) .  We  can  then  simplify 
Eqn.  30  as, 

(I  -  AtD'  +  At  |^(A++  A")  +  At  Bt]AQ  =  -At  Res  (33) 


where  Res  is  the  residual  evaluated  at  time  level  n, 
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,3e  .  3f  3s' 
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Here,  it  is  assumed  that  the  right-hand  side  of  Eqn.  33  (or  Eqn. 

30)  is  evaluated  according  to  Eqn.  20  et  seq.  We  also  note  that 
for  central  differencing  in  rj ,  the  term  incorporating  will 
contain  both  first  and  second  derivatives  and  will  include 
quantities  evaluated  at  j+1,  j,  and  j-1.  The  numerical  solution  of 
Eqn.  33  is  described  in  the  following  subsection. 

Numerical  Solution  Of  The  Thin-Layer  Navier-Stokes  Equations 


The  direct  inversion  of  the  matrix  on  the  left-hand  side  of 
Eqn.  33  contains  unknowns  from  the  five  different  grid-lines  in  the 
£  direction  and  from  three  in  the  i)  direction.  In  particular,  the 
derivative  9(A+AQ)/3£  introduces  unknowns  at  i,  i-1  and  i-2,  while 
9(A_AQ/9£)  introduces  unknowns  at  i,  i+1  and  i+2.  This  results  in 
a  matrix  which  is  quite  costly  to  invert  by  direct  means. 
Consequently,  we  look  to  approximately  factored  procedures  to  solve 
Eqn.  33.  For  this  purpose,  we  express  the  equations  in  fully 
discretized  form  and  factor  using  the  parabolized  Navier-Stokes 
alternating  direction  implicit  (PNS-ADI)  procedure. 

We  can  approximately  factor  the  discretized  version  of  the 
left-hand  side  and  rewrite  Eqn  33  as,  approximately: 


rr  At  . . _  +  A+  .  ,  3 

[r  •  325(4Ai-i‘  Ai-2>  +  3^ 


BTjr'1[r 


At 

2Ax 


<4AI+1-  Ai  +  21AQ  = 


where  the  matrix  T  contains  all  terms  on  the  diagonal, 


r  =  i 


AtD' 


3At 

2Ax 


(36) 


By  performing  the  indicated  multiplication  on  the  left-hand  side  of 
Eqn.  35,  we  see  that  Eqn.  35  is  equal  to  Eqn.  33  except  for  terms 
of  order  At2. 

The  solution  of  Eqn.  35  proceeds  directly.  The  first  operator 
is  tridiagonal  in  T]  with  two  lower  diagonals  corresponding  to  A+i_i 
and  A+i_2*  This  operator  can  be  solved  by  a  sweep  in  the 
streamwise  direction  which  requires  a  block  tridiagonal  matrix 
solution.  The  second  operator  in  brackets  is  lower  (block) 
diagonal  and  can  be  solved  by  marching  from  the  downstream  to  the 
upstream  end  of  the  nozzle.  Thus,  Eqn.  35  represents  a 
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forward-backward  iteration  procedure.  This  iteration  proceeds 
until  AQ  is  driven  to  zero,  or  to  an  acceptable  tolerance.  As  in 
the  perfect  gas  formulation  (see  Appendix  C) ,  this  approximate 
factorization  procedure  is  very  efficient  for  predominantly 
supersonic  flows. 

Solution  Of  The  Euler  Equations 

The  description  above  is  also  applicable  to  the  Euler 
equations,  except  for  several  simplifications.  First  of  all,  the 
viscous  terms  are  omitted  (p  and  k  set  to  zero)  in  Eqns.  13  and  14 
so  the  vector  H'  simplifies  to  only  one  term  and  the  viscous  flux 
S'  vanishes.  Similarly,  the  matrix  in  Eqn.  33  (see  Eqn.  32) 
simplifies  to  B«r  =  B,  and  we  have  a  forward-backward  iteration 
procedure  for  inviscid  flow. 

In  the  special  case  where  the  flow  is  entirely  supersonic,  the 
matrix  A~  becomes  identically  zero,  and  A+  becomes  equal  to  A.  In 
this  case,  the  backward  sweep  degenerates  to  an  identity  operator, 
and  Eqn.  35  reduces  to: 

fr  -  5 h  <4Ai-l  -  At-2>  +  h  B]  A0  =  'at  Res  <37) 

where  T  and  Res  are  simplified  as  noted.  As  in  Appendix  D,  it  is 
recommended  that  each  line  be  time-marched  till  AQ  goes  to  zero 
rather  than  sweeping  the  entire  field.  This  then  results  in  an 
iterative  space-marching  solution  for  supersonic  flows. 

Solution  Of  The  Parabolized  Navier-Stokes  Equations 

The  parabolized  Navier-Stokes  (PNS)  procedure  for  supersonic 
flow  is  based  upon  removing  the  terms  in  the  equations  that  allow 
upstream  propagation  of  information.  This  involves  two  steps.  The 
first  step  involves  removing  viscous  diffusion  in  the  streamwise 
direction  which  has  already  been  done  in  the  TLNS  equations.  The 
second  step  involves  omitting  the  upstream  propagation  of 
information  through  the  subsonic  portion  of  the  boundary  layer. 

These  approximations  make  the  equations  parabolic  and  allow  space 
marching  procedures  to  be  used. 

The  present  parabolized  procedure  is  obtained  by  setting  the 
matrix  A~  to  zero  in  the  TLNS  formulation  of  Eqn.  35.  By  omitting 
A",  the  back  sweep  of  Eqn.  35  simplifies  to  the  identity  matrix  and 
the  remaining  algorithm  corresponds  exactly  to  the  forward  sweep  of 
the  full  TLNS  procedure.  In  omitting  A“,  we  drop  all  upstream 
running  characteristics  and  thus  obtain  a  well-posed  space  marching 
procedure.  The  resulting  numerical  algorithm  becomes  identical  to 
the  inviscid  procedure  for  the  supersonic  Euler  equations  except 
that  the  viscous  terms  are  retained.  (Recall  the  A-  operator  is 
identically  zero  in  supersonic  flow.)  Setting  A"  to  zero  and 
eliminating  the  upstream  operator  is  completely  analogous  to 
multiplying  the  pressure  gradient  by  a  value  of  u)  less  than  one  and 
ignoring  the  (l-w)3p/3£  term  as  is  done  in  standard  PNS  procedures5. 
The  advantage  of  the  present  procedure  is  that  it  is  more  firmly 
founded  from  a  physical  viewpoint. 


(38) 


Mathematically,  the  PNS  algorithm  is  expressed  as, 

‘r  -  5S5  <4At-l  -  AI-2>  +  h  -  At  Res 

where  T  is  as  defined  in  Eqn.  36  except  that  A~  is  zero,  and  B<j> 
still  retains  its  full  viscous  form  of  Eqn.  32.  As  for  the 
supersonic  inviscid  flow  calculation,  iterations  at  each  line  are 
converged  before  moving  downstream  so  that  the  converged  solution 
is  obtained  in  a  single  sweep. 

Euler  Equation  Boundary  Conditions 

Boundary  conditions  are  specified  in  exactly  the  same  method 
as  described  in  Appendix  C.  We  begin  with  boundary  conditions  for 
the  Euler  equations.  For  the  Euler  equations,  we  rely  on  the 
Method  of  Characteristics  to  ensure  we  enforce  both  the  correct 
number  of  conditions  and  types  of  conditions  that  are  physically 
and  mathematically  acceptable.  In  general,  this  means  that  the 
four  unknowns  (corresponding  to  the  four  equations  in 
two-dimensional  flow)  are  obtained  from  a  combination  of  boundary 
conditions  and  specific  subsets  of  the  equations  of  motion. 

At  a  subsonic  inflow  boundary,  we  specify  three  quantities  as 
boundary  conditions  and  obtain  the  last  from  the  equations  of 
motion.  For  supersonic  inflow,  we  specify  all  four  quantities. 
Since  we  are  primarily  interested  in  supersonic  nozzle  flows  which 
are  started  downstream  of  the  throat  of  a  converging-diverging 
nozzle,  this  latter  choice  is  the  one  of  primary  interest  here. 
Consequently,  we  give  no  details  as  to  the  MOC  procedure  for 
subsonic  inflow.  These  can,  however,  be  obtained  from  an  extension 
of  the  MOC  procedures  for  wall  boundary  conditions,  or  from 
Appendix  C.  In  the  case  of  supersonic  inflow,  the  necessary  two 
lines  of  information  are  obtained  from  the  transonic  solution  (see 
Appendix  C) . 

At  the  outflow  boundary,  the  solution  of  Euler  equations  for  a 
diverging  nozzle  is  wholly  supersonic.  Therefore,  all  the 
information  at  the  downstream  end  is  determined  from  the  flowfield 
itself.  No  downstream  boundary  condition  is  needed  or  allowed.  On 
the  axis  of  symmetry,  y  =  0,  therefore, 


Q  =  0 


(39) 


because  y  is  contained  inside  Q  (see  Eqn.  1) .  After  Q  is 
determined,  the  primitive  variables  are  determined  by  a  second 
order  extrapolation.  If  f  denotes  a  primitive  variable,  then  on 
the  axis  9f/9y  =  0  and, 


f 


5yfi-~l  ~  _ ( 4  f  1  ,  2  f  i  ,  3 ) 

Sy  +  S/2ny 


(40) 


Since  v  =  0  on  the  axis,  there  are  only  three  primitive  variables 
to  be  extrapolated.  The  real  gas  variables  are  f  =  p,  u,  eQ. 

For  inviscid  flows,  slip  wall  boundary  conditions  are  enforced 
by  means  of  the  Methods  of  Characteristics  and  auxiliary 
information  corresponding  to  in-running  characteristics. 
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From  an  eigenvalue  analysis  of  the  Jacobian,  B  =  3f/3t),  one 

gets, 

Ab  =  diag  (V,  V,  V+C,  V-C) 

where  V  is  the  contravariant  velocity  component,  V  =  twu  +  r|yV,  and 
C  is  the  transformed  speed  of  sound,  C  =  c(i]x2  +  Tly2)1'2. 

On  a  wall,  the  first  three  eigenvalues  are  positive  or  zero 
and  the  fourth  is  negative.  Accordingly,  three  quantities  on  the 
boundary  are  enforced  by  outgoing  characteristic  relations  and  the 
fourth  is  obtained  from  the  specified  boundary  condition,  V  =  0 
which  forces  the  flow  to  be  tangent  to  the  wall. 

The  three  characteristic  equations  are  derived  from  the 
discretized  Euler  equations  (Eqns.  30)  by  premultiplying  both  the 
left-  and  right-hand  sides  by  the  product  of  a  selection  matrix,  L, 
and  an  eigenvector  matrix,  Tg“lp-1,  where  the  matrix  L,  L  =  diag 
(1,  1,  1,  0) ,  acts  as  a  selection  operator  that  selects  only  the 
information  (equations)  corresponding  to  the  outgoing 
characteristics.  The  modal  matrix,  Tg”lp~l  is  given  in  Appendix  A. 

The  single  boundary  condition  which  combines  with  the  three 
characteristic  relations  given  above  can  be  expressed  in  vector 
form  as, 


fl(Q)  =  [0,  0,  0,  V]T  =  0. 

In  order  to  express  this  boundary  condition  in  delta  form,  we 
expand  Q  in  a  Taylor's  series, 

nn+l  _  nn  +  3Q  iQ  _  0  ,41) 

where  30/3q  is  a  Jacobian  matrix. 

By  combining  the  one  relation  in  Eqn.  41  with  the  three 
conditions  in  Eqn.  40,  we  obtain  four  coupled  equations  which 
determine  the  four  unknowns  on  the  nozzle  wall. 

Viscous  Boundary  Conditions  For  The  TLNS  Equations 

For  the  thin-layer  Navier-Stokes  equations,  two  boundary 
conditions,  namely  those  at  the  upstream  end  and  those  along  the 
axis,  are  identical  to  those  applied  to  the  Euler  equations. 
Additional  detail  is  not  given  here.  Both  the  outflow  boundary 
conditions  and  the  solid  wall  boundary  conditions  are  different. 

The  outflow  conditions  are  still  treated  in  an  inviscid  manner  by 
the  Method  of  Characteristics,  but  because  of  the  nozzle  wall 
boundary  layer,  portions  of  the  flow  at  the  nozzle  exit  are 
subsonic,  so  both  subsonic  and  supersonic  conditions  are  needed 
there.  At  the  wall,  the  viscous  no  slip  conditions  are  specified. 
Details  of  these  implementations  are  given  here. 

The  flow  at  the  downstream  boundary  consists  of  two  parts  —  a 
major  supersonic  part  extending  from  the  axis  to  the  vicinity  of 
the  wall,  and  a  small  subsonic  part  adjacent  to  the  wall 
corresponding  to  the  subsonic  part  of  the  boundary  layer. 
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Quantities  on  these  two  parts  of  the  boundary  line  are  determined 
differently  because  of  the  different  nature  of  supersonic  and 
subsonic  flows.  The  fourth  eigenvalue  of  the  Jacobian,  A  =  9e/3q, 
changes  sign  from  positive  to  negative  as  the  flow  changes  from 
supersonic  to  subsonic. 

From  an  eigenvalue  analysis  of  the  Jacobian  A,  one  gets, 

QA  =  diag  (U,  U,  U+C,  U-C) 

where  U  is  the  contravariant  velocity  component,  U  -  £xu  +  £yv,  and 
C  is  the  transformed  speed  of  sound,  C  =  R(£x2  +  In  the 

supersonic  region,  all  eigenvalues  are  positive.  Therefore,  the 
boundary  conditions  are  determined  by  the  outgoing  flow  which  is 
solved  by  the  forward  sweep  step.  In  the  subsonic  region,  the 
first  three  eigenvalues  are  positive  and  the  fourth  is  negative. 
This  means  that  three  of  the  four  unknowns  are  determined  by  the 
outgoing  characteristics  while  the  fourth  must  be  imposed  as  a 
boundary  condition. 

The  three  characteristic  equations  are  derived  from  the 
discretized  TLNS  equations  (Eqn.  33)  upon  premultiplying  by 
LTa-I-P-I,  in  the  same  manner  as  described  for  the  inviscid  wall 
boundary  condition.  The  difference  is  that  here  the  selection 
matrix  is  given  by  L  =  diag  (1,  1,  1,  0)  so  that  only  the 
information  (equations)  corresponding  to  the  outgoinq 
characteristics  is  retained.  A  second  difference  is  that  the  modal 
matrix,  TA~lp~l  (given  in  Appendix  A)  is  the  one  corresponding  to 
the  Jacobian  A. 

The  boundary  condition  that  is  used  to  augment  these  three 
Method  of  Characteristic  relations  is  a  specification  of  the  nozzle 
back  pressure,  pb,  over  those  portions  of  the  flow  that  are 
subsonic.  In  vector  notation,  this  becomes, 

fi(Q)  =  [0,  0,  0,  p  -  pb]T  =  0  (42) 

Then,  taking  the  Jacobian  of  d  with  respect  to  Q,  we  obtain  as 
before, 

^  AQ  =  pb-  p  (43) 

The  right-hand  side  of  Eqn.  43  can  either  be  left  as  is  to  prevent 
small  drifts  of  the  pressure  at  the  exit  plane,  or  can  be  set  to 
zero  as  Eqn.  42  suggests. 

In  addition  to  this  mathematically  and  physically  correct 
treatment  of  the  subsonic  portion  of  the  outflow,  we  also  use  a 
second  option  of  extrapolating  the  subsonic  pressure  from  inside 
the  flowfield.  This  second  procedure  is  not  correct 
mathematically,  but  it  is  frequently  imposed  in  problems  of  this 
type.  Its  success  rests  with  the  smallness  of  the  subsonic  region. 
Choice  of  the  extrapolation  condition  eliminates  the  possibility  of 
determining  the  effect  of  perturbations  in  the  back  pressure  on 
nozzle  performance,  an  issue  that  is  sometimes  of  considerable 
interest. 


To  apply  the  extrapolated  pressure  condition,  we  use  the  three 
characteristic  relations  along  with  a  first  order  extrapolation  of 
pressure,  namely, 


P 


n+1 


n+1 


(44) 


where  I  is  the  last  grid  line.  In  vector  form,  this  is  expressed 
as. 


f)j  -  Qj-l  =  0  (45) 

where , 

»i  =  (0,  0,  0,  Pij)T  (46) 

Finally,  we  note  that  in  using  the  specified  boundary 
condition,  pg,  only  minor  deviations  from  the  ideally  expanded 
nozzle  can  be  accommodated.  A  typical  procedure  is  to  first 
determine  the  "ideally  expanded"  back  pressure  by  the  extrapolation 
condition.  Then  other  back  pressures  in  the  neighborhood  of  this 
condition  can  be  specified.  Generally  speaking,  the  back  pressure 
can  only  be  varied  by  a  factor  of  two  or  three  from  the  value 
determined  by  the  extrapolation  procedure.  This  variation  is, 
however,  enough  to  show  a  decided  thinning  of  the  subsonic  portion 
of  the  boundary  layer  when  the  back  pressure  is  lowered,  as  well  as 
a  thickening  of  the  subsonic  region  and  the  establishment  of  a 
small  recirculating  region  when  the  back  pressure  is  increased. 

When  re-entry  flow  is  encountered,  inflow  conditions  must  be  added 
for  those  points  on  the  exit  boundary  that  experience  inflow.  In 
general,  the  procedure  fails  when  the  inflow  region  becomes  large 
enough  that  the  point  of  demarcation  between  inflow  and  outflow 
begins  to  jump  back  and  forth  during  the  iteration.  To  go  beyond 
these  limits,  it  is  necessary  to  expand  the  flow  domain  to  extend 
outside  the  nozzle.  Nevertheless,  the  present  procedure  can  be 
very  effective  in  understanding  the  physics  of  non-ideally 
expanded,  viscous  supersonic  flows  and  in  assessing  performance 
shifts  when  the  nozzle  back  pressure  is  nearly,  but  not  exactly, 
matched  to  the  flow. 

For  the  viscous  wall  boundary  condition,  we  specify  three 
boundary  conditions.  The  no-slip  boundary  condition  implies  that 
both  components  of  the  contravariant  velocity  are  zero,  u  =  V  =  0. 
The  third  boundary  condition  is  one  on  the  energy  flux  at  the  wall. 
This  can  either  take  the  form  of  a  specified  wail  temperature, 

(Tw  =  Tw(£) ,  or  a  specified  heat  flux,  8T/3y  =  f(£),  where  f  is  an 
arbitrary  function. 

The  fourth  condition  for  the  viscous  wall  boundary  condition 
comes  by  solving  the  radial  momentum  equation  for  the  normal 
pressure  gradient  at  the  wall,  and  enforcing  this  momentum  equation 
in  conjunction  with  the  three  boundary  conditions.  The  normal 
gradient  of  the  pressure  in  a  non-orthoqonal  coordinate  system  is 
given  by, 


n*Vp  =  (T)xe;x  -  T)y£y)  P£  +  (!)x2  +  T) y 2  ;  pi)  (47) 
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where  n  is  the  unit  normal  to  the  wall.  An  analogous  relation  is 
used  for  specifying  the  temperature  gradient  when  the  heat  flux 
boundary  condition  is  imposed.  In  differencing  Eqn.  47,  second 
order  accurate  one-sided  differences  are  used  in  r]  and  central 
differences  are  used  in  £. 

Boundary  Conditions  For  The  PNS  Procedure 

Boundary  conditions  for  the  PNS  procedure  are  identical  to 
those  for  the  TLNS  equations  except  that  no  downstream  conditions 
are  allowed  and  no  downstream  information  is  allowed  in  computing 
normal  derivatives  on  the  walls.  When  normal  derivatives  of 
pressure  and  temperature  are  obtained,  similar  procedures  are  used 
except  the  £  derivatives  are  evaluated  by  a  second-order  backward 
difference. 

Computing  The  Rate  Of  Mass  Flow  And  Thrust 

For  an  axisymmetric  nozzle,  the  conservation  of  mass  flow  is 
verified  by  integrating  along  arbitrary  cross-sections  at  various 
streamwise  positions.  With  no  mass  addition,  the  rate  of  mass  flow 
must  be  equal  at  each  of  these  sections.  The  nozzle  thrust  is 
determined  by  integrating  across  the  cross-section  at  the  exit 
plane . 

For  simplicity,  cross-sections  along  constant  values  of  £  are 
chosen  to  verify  the  conservation  of  mass  flow  and  to  compute  the 
thrust . 

The  rate  of  mass  flow,  dm,  through  an  elemental  surface,  ds, 
is  given  by, 

dm  =  pU£ds  (48) 

where, 

U  -  - y -  =  - £xH_+ — LyY —  (49) 

E  <£x2  ♦  EyV/2  (Cx2  ♦  ty2>1/2 

is  the  component  of  the  velocity  vector  normal  to  the  surface, 

£  =  £(x,y),  and  ds  =  2nydy,  is  the  elemental  area  at  a  distance  y 

from  the  axis. 

In  Cartesian  notation,  dr)  =  ((dx)2  +  (dy)2)l/2f  and  the  rate 
of  mass  flow  is  given  by, 

dm  =  2 n  (  puydy  -  pvydx)  (50) 

The  integration  of  dm  is  performed  numerically  between  j=l  and  j=JL 
for  each  i  using  Simpson's  rule. 

In  the  real  gas  TLNS  equations,  the  fluxes  are  defined  at  the 
midpoint  Ei+i/2j  (Eqn.  25).  Then,  for  second  order  accuracy  (the 
index  j  is  omitted  for  simplicity), 


(PUf;)^  =  5  f  CPUe)i  +  (PUf)i+1l 

-  |  1*1  (P«t)|  i+k~  ^(PVU  +  a<PU)I+3/2)  <51> 

The  differential  thrust,  dt,  in  the  axial  direction,  is  given 

by, 

dip  =  dmu  +  (P  -  Pfc)  2nydy  (52) 

where  Pfc>  is  the  prescribed  pressure  at  the  exit. 
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USERS  MANUAL  FOR  THE  CODES 


Introduction 

The  following  discussion  is  intended  to  be  an  overview  that 
explains  briefly  how  to  run  the  codes.  A  total  of  four  different 
codes  are  described  that  solve  the  perfect  or  real  gas  flow  for  the 
Euler  or  thin-layer  Navier-Stokes  equations  in  an  axisymmetric 
nozzle.  In  particular,  codes  AXI2DSA,  AXI2DSE  and  AXI2DSF  solve 
the  Euler  equations,  while  the  codes  DDADIPBC,  DDADIPBE  and 
DDADIPBG  solve  the  TLNS  equations. 

The  perfect  gas  flow  is  solved  by  the  codes  AXI2DSA,  AXI2DSE, 
DDADIPBC,  and  DDADIPBE.  The  real  gas  flow  is  solved  by  the  codes 
AXI2DSG  and  DDADIPBG.  The  real  gas  codes  are  demonstrated  for 
equilibrium  air  flow,  however,  by  changing  the  input  data  and  the 
relevant  data  subroutines,  the  user  is  able  to  solve  any  other 
equilibrium  real  gas  flow.  These  modifications  are  explained  later 
in  this  section. 

The  Euler  solver,  AXI2DSA,  and  the  TLNS  solver,  DDADIPBC,  are 
designed  to  solve  perfect  gas  flows.  In  these  codes,  the  molecular 
weight,  Mw,  is  constant,  and  the  internal  energy,  e,  is  given  by, 


e  =  cvT  (53) 

where  cV/  the  specific  heat  coefficient  in  constant  volume,  is 
constant,  and  T  is  the  temperature.  In  the  TLNS  solvers  the 
laminar  viscosity  coefficient,  p,  may  be  defined  by  one  of  the 
following  three  options: 


a) 

b) 


c) 


constant  value,  p  =  pQ 
Sutherland  law, 


u  =  u  +  196  T_ 
M  Mo  T  +  196  'T0' 


1.5 


Power  law, 


u> 


v  =  »*o  (T^> 


(54a) 


(54b) 


( 54c) 


where  T0  and  pQ  are  reference  values  and  w  is  a  constant  ranging 
between  0.5  £  w  ~  1.0. 

The  Euler  code  AXI2DSC  and  the  TLNS  code  DDADIPBG  are  the  most 
generalized  solvers.  In  these  codes,  the  molecular  weight,  Mw,  and 
the  internal  energy,  e,  are  defined  as  functions  of  two  variables. 

Mw  =  Mw  (p,T)  (55) 

and 


e  =  e(p,T),  or  T  =  T(p,e) 


(56) 
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In  literature  these  relations  appear  as  tables6'7.  The  viscosity 
and  thermal  conductivity  coefficients,  y  and  k,  are  also  specified 
as  functions  of  two  variables, 

P  =  P(P#  e) 
k  =  k(p,  e) 

in  the  real  gas  TLNS  code.  Simplified  curve  fits  for  y  and  k  for 
air  were  suggested  by  Srinivasan  et  al.8. 

In  order  to  run  the  codes,  the  user  must  furnish  the 
transformed  coordinates  for  the  flowfield  and  define  the  gas 
properties.  The  real  gas  properties,  excluding  the  viscosity  and 
thermal  conductivity,  are  read  into  the  Euler  and  TLNS  codes.  The 
viscosity  and  thermal  conductivity  data  are  generated  by  the  TLNS 
codes  during  the  computation  process. 

In  the  Euler  codes,  AXI2DSA  and  AXI2DSG,  the  flowfield  initial 
guess  at  each  i  and  i  =  1,2,3...,  JL  is  obtained  from  the  previous 
converged  solution  at  i-1.  Since  the  Euler  solution  converges 
relatively  rapidly,  it  is  usually  obtained  in  a  single  run. 
Consequently,  the  constants  that  monitor  the  number  of  iterations 
and  the  restart  for  the  next  run,  although  defined  in  namelist 
INPUT,  are  not  in  use. 

In  the  TLNS  codes,  DDADIPBC  and  DDADIPBG,  the  flowfield 
solution  is  obtained  by  alternate  forward  and  backward  sweeps.  The 
flowfield  initial  guess  is  evaluated  by  two  ways.  In  DDADIPBC,  the 
initial  guess  is  obtained  by  a  PNS  single-sweep  procedure  which  is 
an  integrated  part  of  the  code.  In  DDADIPBG,  the  flowfield  initial 
guess  is  derived  from  the  isentropic  relations  throughout  the  flow 
domain.  The  above-described  PNS  initial  condition  represents  a 
cost  savings  of  two  or  more  and  its  implementation  in  this  code  is 
recommended . 

The  iterative  solution  of  the  TLNS  codes  requires  a  large 
amount  of  computer  time.  Therefore,  the  user  should  be  cautious  in 
defining  the  constants  that  monitoring  the  number  of  iterations  and 
continuation  of  the  next  runs . 

Table  2  presents  a  summary  of  the  codes'  executing  data,  which 
includes  the  number  of  mesh  points,  the  recommended  number  of 
iterations,  and  the  computing  time  per  iteration  per  mesh  point  on 
the  NOS-BE  system.  In  this  table,  the  number  of  mesh  points 
denotes  the  maximum  amount  allowed  by  the  core-memory  of  the  CYBER 
180/840,  NOS-BE  system.  As  indicated  before,  these  stringent 
limitations  are  removed  when  running  under  NOS-VE.  The  converged 
solution  is  obtained  when  the  average  relative  error,  A,  defined 
as, 


A 


1 

JL*  IL 


JL 

IL 

E 

E 

j=i 

i=l 

Mu4 

Qilj 


(57) 


is  less  than  a  specified  value.  In  the  present  runs,  this  was  set 
at  A  =  lO-!®  for  the  Euler  codes  and  A  =  10-6  for  the  TLNS  codes. 

The  units  used  in  the  codes  are  the  International  System  of 
Units  (SI  units).  These  SI  units  are  as  given  in  Table  3.  The 
universal  gas  constant  that  is  consistent  with  these  units  is 
R  =  8314.3  J/ (kgmole) «K. 
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In  the  energy  equation,  the  units  of  the  internal  energy,  e, 
are  given  by  (m/sec)2. 

Additional  Information  Concerning  Real  Gas  Flows 

Calculations  of  real  gas  functions  require  the  molecular 
weight,  the  temperature,  the  viscosity  and  thermal  conductivity  and 
several  derivatives  of  these  quantities  to  be  known  as  functions  of 
two  thermodynamic  variables.  As  primary  dependent  variables,  we 
choose  the  properties  density,  p,  and  internal  energy,  e.  A  number 
of  FUNCTION  subroutines  have  been  written  to  obtain  these 
quantities  from  tabular  input  data.  These  functions  are  described 
in  Table  4.  In  addition  to  these  Function  Subprograms,  subroutines 
UGAS3  and  UGAS4  are  called  by  Functions  FZMU  and  FZK. 

The  Perfect  Gas  Euler  Code,  AXI2DSA 

The  code  AXI2DSA  solves  the  Euler  equations  for  a  perfect  gas 
flow.  The  present  version  has  all  dimensional  arrays  set  to  (70  x 

44),  i.e.,  70  points  in  the  £  direction  and  44  points  in  the 

direction.  These  are  the  maximum  grid  points  allowed  by  the  CYBER 

NOS-BE  memory.  These  limitations  can  easily  be  reset  by  a 

recompilation  of  the  code.  Only  the  DIMENSION  statements  need  be 
changed,  all  other  FORTRAN  statements  will  execute  properly  if 
these  numbers  are  increased. 


Input  Data  Description 


1.  Namelist  INPUT.  Namelist  INPUT  specifies  the  flow 
constants.  These  constants  are  read  from  file  TAPE1  by  subroutine 
INITIA.  The  quantities  in  this  NAMELIST  are  as  follows: 

IL  -  Number  of  mesh  points  in  the  £  direction  (Integer) . 

JL  -  Number  of  mesh  points  in  the  T]  direction  (Integer)  . 

ITR1  -  Not  used;  set  to  1. 

ITRN  -  Not  used;  set  to  1. 

ITRS  -  =  ITRN  -  ITR1  +  1  (Integer) . 

IS VP  -  Not  used;  set  to  0. 

CFL  -  Courant-Friedrichs-Levy  number  (typical  ranges 
100-500) . 

ITIME  -  Time  step  control  (Integer)  : 

ITIME  =  0  for  constant  At; 

ITIME  =  1  for  constant  CFL  (normal  value) . 

THETA  -  Crank-Nicholson  constant,  always  set  to  1.0. 

RM1  -  The  Mach  number  at  the  throat,  RM1  >  1.0. 

RM2  -  The  exit  Mach  number;  approximate  value,  only  used  for 
set-up. 

AIN  -  Not  used. 

AEX  -  Not  used. 

RL  -  Not  used. 

PO  -  Upstream  stagnation  pressure,  n/m2  (uniform  across  the 

stream) . 

TO  -  Upstream  stagnation  temperature,  K  (uniform  across  the 
stream) . 

CP  -  Specific  heat  in  a  constant  pressure  (constant) . 
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cv 

GAMMA  - 
GM1 

INER  - 

IB1 

IBL 

NORD  - 


OMEGAX- 

OMEGAY- 

IREAD  - 
IWRT  - 

IRVN  - 


Specific  heat  in  a  constant  volume  (constant) . 

Ratio  of  specific  heats,  CP/CV. 

(GAMMA-1) . 

Maximum  Number  of  iterations  for  a  converged  solution 
(Integer) . 

Grid  line  in  the  £  direction  at  which  the  iterative 
process  starts.  Usually  IBl  =  2  (Integer). 

Grid  line  in  the  u  direction  at  which  the  iterative 
process  ends.  Usually  IBL  =  IL  (Integer) . 

Control  variable  for  specifying  the  order  of  accuracy 
in  the  £  direction  (Integer) ; 

NORD  =  0  -  first  order  accurate 

NORD  =  1  -  second  order  accurate  (normal  value) . 
Artificial  dissipation  in  the  £  direction.  (Real 
variable,  generally  less  than  0.5). 

Artificial  dissipation  in  t|  direction.  (Real  variable 
generally  less  than  0.5). 

Not  used;  set  to  0. 

Output  print  control  (Integer) : 

IWRT  =  o  -  Brief  printout 
IWRT  =  1  -  Full  printout 
Not  used;  set  to  1. 


2.  Mesh  Points  Coordinates.  The  flowfield  mesh  points 
coordinates,  X(I,J),  Y(C,J)  for  1=1 , 2 , 3 , . . . , IL,  J=1 , 2 , 3 , . . . , JL,  are 
read  from  file  TAPE 3  by  subroutine  INITIA.  In  this  input  data  the 
inflow  boundary  conditions  are  defined  on  1=1,  J=1 , 2 , 3 , . . . , JL,  the 
outflow  on  I=IL,  J=l, 2, 3, . . . ,TL,  the  axis  of  symmetry  on  J=l, 

1=1 , 2 , 3 , . . . , IL,  and  the  wall  on  J=JL,  1=1 , 2 , 3 , . . . , IL.  The  reading 
format  is, 

READ  (3,401) ( (X(I,J) ,  Y(I,J),  1=1, IL),  J=1,JL) 

501  FORMAT  (E17.9,  4E16.9) 


Output  Data  Description 

1.  Printed  Results.  For  the  short  form  printout,  IWRT  is 
set  to  0  in  namelist  INPUT.  With  this  printout  choice,  the 
convergence  rate  for  each  equation  from  subroutine  SUPPLY  (ENTRY 
CHECK) ,  and  the  downstream  rate  of  mass  flow  from  subroutine  SUPPLY 
(ENTRY  MASS)  is  printed  at  each  iteration. 

When  IWRT=1  a  detailed  printout  of  the  flowfield  results  is 
printed  by  subroutine  SUPPLY  (ENTRY  OUTPUT) .  For  each  mesh  point 
(I,J)  it  includes  the  coordinates,  (X,Y),  the  velocity  components 
(u,v) ,  the  pressure,  density,  temperature,  stagnation  internal 
energy,  entropy,  and  the  Mach  number. 

The  printout  is  stored  on  TAPE6  and  printed  automatically  at 
the  end  of  the  run. 


2.  Conservative  Variables.  For  further  evaluation  such  as 
plotting,  the  iteration  constants,  ITR1,  ITRN,  the  local  time 
increment,  DELTAU(I,J)  (for  Atfij),  and  the  conservative  variables, 
RHO ( I , J ) ,  RHOU ( I , J ) ,  RHOV ( I ,  J )  and  EO(I,J)  (for  Piij,  (pu)iij, 
(Pv)iijf  and  (eo)iij'  respectively),  are  written  on  TAPES  by 
subroutine  SUPPLY  (ENTRY  OUTPUT),  as  follows, 
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WRITE  (8,504)  ITRl ,  ITRN 
504  FORMAT  (IG,  IG) 

WRITE  (8,503)  ( (DELTAV ( I , J) ,  1=1, IL),  J=1,JL) 

503  FORMAT  (E17.9,  4E16.9) 

WRITE  (8,503)  ((RHO(I,J),  RHOU(I,J),  RHOV(I,J), 

E (I , J) ,  1=1, IL) ,  J=1 , JL) 


The  Real  Gas  Euler  Code,  AXI2DSG 

This  code  solves  the  Euler  equations  for  a  real  gas  flow.  The 
molecular  weight,  Mw,  and  the  internal  energy,  e,  are  functions  of 
two  variables,  p  and  T.  For  air,  these  relations  are  implemented 
by  the  code  INPLA  and  then  are  read  as  an  input  by  the  mail  code 
AXI2DSG. 

In  the  present  version  of  AXI2DSG,  the  number  of  mesh  points 
in  the  £  and  t]  directions  are  set  at  70  and  44,  respectively. 


Input  Data  Description 

1.  Namelist  INPUT.  This  namelist  is  identical  to  that  for 
the  perfect  gas  code  described  in  the  previous  section  except  that 
the  specific  heat  inputs  CP  and  CV  along  with  the  specific  heat 
ratios  GAMMA  and  GMI  are  omitted.  In  their  place,  the  following 
gas  constant  inputs  are  added: 

RG  -  Universal  Gas  Constant  (8314.3). 

AMWO  -  Gas  molecular  weight  at  inlet  stagnation  conditions 
(uniform  across  stream) . 

GAMMAO-Ratio  of  specific  heats  at  inlet  stagnation 
conditions  (uniform  across  stream) . 


2.  Namelist  DINPL.  This  namelist  reads  in  the  real  gas 
properties  for  the  gas  of  interest.  The  data  are  read  from  TAPE 9 
by  subroutine  INITIA. 

XI (I)  -  Array  of  derivatives  for  molecular  weight. 

Y1(J)  -  Array  of  temperatures  for  molecular  weight. 

F1(I,J)-  Array  of  molecular  weights  at  (XI (I),  Yl(J)). 


X2 ( I )  -  Density  array  for  molecular  weight  derivative. 

Y2(J)  -  Temperature  array  for  molecular  weight  derivative. 

F2 ( I , J)  -  Array  of  molecular  weight  derivatives,  (9Mw/9p)i>  at 
(X2(I),  Y2 ( J) ) . 


X3 (I)  -  Density  array  for  9mw/9t. 

Y3(J)  -  Temperature  array  for  (9mw/9t) . 

F3 (I , J) -  Array  of  the  derivatives  (9mw/9t)p  at  (X3(I), 
Y3 (J) )  . 


XY(F)  -  Density  array  for  internal  energy. 

Y4(J)  -  Temperature  array  for  internal  energy. 

F4 ( I , J) -  Array  of  internal  energy  values,  e,  at  (XY(I), 
Y4 ( J) )  . 


X5 (I)  -  Density  array  for  temperature. 

Y5(J)  -  Internal  energy  array  for  temperature. 
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F5 (I, J) -  Array  of  temperatures  at  (X5(I),  Y5(J). 

X6 (I)  -  Density  array  for  3T/3p 

Y6(J)  -  Internal  energy  array  for  9t/3o 

F6 (I , J) -  Array  of  derivatives,  (3T/3p)e  at  (X6(I),  Y6(J)). 

X7 ( I )  -  Density  array  for  9T/9e. 

Y7(J)  -  Internal  energy  array  for  3T/3e. 

F7 (I , J) -  Array  of  derivatives,  (3T/9e)p  at  (X7(I),  Y7(l)). 

3.  Mesh  Points  Coordinates.  The  transformed  coordinates 
are  read  in  the  same  manner  as  for  the  perfect  gas  code.  See  the 
section  on  The  Perfect  Gas  Euler  Code,  AXI2DSA. 

Output  Data  Description 

1.  Printed  Results.  For  the  shortform  printout,  IWRT=0, 
the  namelists  INPUT  and  DINPL,  the  iterated  values  of  the  molecular 
weight,  Mw,  and  the  effective  Y  at  the  throat  are  printed  from 
subroutine  INITIA.  In  addition,  the  convergence  rate  for  each 
equation  is  printed  from  subroutine  SUPPLY  (ENTRY  CHECK) ,  and  the 
downstream  rate  of  mass  flow  is  printed  from  subroutine  SUPPLY 
(ENTRY  MASS) . 

The  print  for  IWRT=1  gives  all  flowfield  quantities  listed  in 
Output  Data  Description  for  the  perfect  gas  Euler  code.  In 
addition,  it  prints  the  molecular  weight  and  the  entropy.  The 
printout  is  again  stored  on  TAPE 6  and  printed  automatically  at  the 
end  of  the  run. 


2.  Conservative  Variables.  Plotting  information  is  also 
given  on  TAPES  in  the  same  format  and  for  the  same  variables  as 
listed  in  Output  Data  Description  for  the  perfect  gas  Euler  code. 

The  Perfect  Gas  TLNS  Code,  DDADIPBC 


This  code  solves  the  TLNS  equations  for  a  perfect  gas  flow. 

The  laminar  viscosity,  p,  is  determined  by  the  relations  in  Eqns. 
54a-c  and  the  turbulent  viscosity,  pf,  is  determined  by  the 
Baldwin-Lomax  model  of  turbulence^.  The  molecular  thermal 
conductivity,  k,  and  the  turbulent  thermal  conductivity,  k<p,  are 
expressed  in  terms  of  the  laminar  turbulent  viscosities  by  means  of 
the  Prandtl  numbers,  Pr,  and  Pr>p,  respectively,  as  defined  in  Eqn. 
13. 


In  the  perfect  gas  viscous  code,  the  maximum  number  of  mesh 
points  as  limited  by  the  NOS-BE  system  are  (70  x  44) ,  the  same  as 
for  the  perfect  gas  Euler  code.  Because  of  the  additional  run  time 
required  (see  Table  2) ,  the  viscous  code  is  designed  to  be  stopped 
and  restarted. 


Input  Data  Description 

1.  Namelist  INPUT.  As  for  the  inviscid  codes,  this 
namelist  specifies  the  flow  constants.  The  namelist  is  read  from 
TAPEl  by  subroutine  INITIA.  Most  of  the  variables  are  analogous  to 
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those  described  in  Input  Data  Description  for  the  perfect  gas  Euler 
code.  These  have  the  same  definitions  as  given  there.  Additional 
quantities  for  specifying  the  viscous  flow  are: 

NBEG  -  The  initial  iteration  number.  Generally  set  to  1. 

NEDN  -  Total  number  of  iterations  to  be  performed.  If 
NEND=1,  the  PNS  solution  is  obtained. 

NITER  -  Number  of  iterations  in  the  run  (Integer) . 

CFL1  -  CFL  number  for  the  PNS  solution. 

IVISC  -  An  integer  specifying  the  nature  of  the  flow: 

IVISC=0  inviscid  flow  (with  IVISC=0,  this  code 
should  duplicate  the  results  of  AXI2DSA) ; 

IVISC=1  viscous  flow. 

IWALL  -  An  integer  specifying  the  wall  boundary  condition; 
IWALL=0  an  adiabatic  wall, 

IWALL=1  specified  wall  temperature. 

PRN  -  Laminar  Prandtl  number  (constant) . 

PRNT  -  Turbulent  Prandtl  number  (constant) . 

TREF  -  The  reference  temperature  for  the  viscosity  in  Eqn. 
54,  K. 

ZMUO  -  Reference  viscosity  at  T=TREF  (SI  units) . 

OMEGA  -  Exponent  for  the  power  law  form  of  the  viscosity, 
Eqn.  54c. 

TWALL  -  Value  of  the  wall  temperature  when  IWALL=1 
(TWALL=0 . 0  when  IWALL=0) ,  K. 

PB  -  The  back  pressure  at  the  nozzle  exit.  When  PB=0, 
the  subsonic  pressure  at  the  exit  is  extrapolated 
from  interior. 

IRVN  -  Control  variable  for  restart  (Integer) : 

IRVN=0,  for  the  first  run  only, 

IRVN=A,  for  restart  from  previous  runs. 

2.  Mesh  Points  Coordinates.  These  are  read  in  the  same 
format  as  for  the  Euler  code,  Input  Data  Description. 

3.  Convergence  Rate,  Previous  Run.  (Ignored  if  not  a 
restart  case.)  For  the  case  of  restart  files,  additional 
(formatted)  input  is  placed  on  TAPE1  to  specify  the  iteration 
number  (NDUM+1)  at  which  the  restart  run  begins.  This  input  is 
also  used  to  create  a  continuous  file  of  the  convergence  of  the 
TLNS  equations.  This  input  is  placed  on  TAPE1  behind  the  namelist 
INPUT  data.  It  is  also  read  by  subroutine  INITIA.  The  reading 
format  is, 

70  READ  (9,  502,  END=65)  NDUM,  (SS(K),  K=1,Y) 

502  FORMAT  (15,  3X,  4(1X,  E14.7) 


GO  TO  70 
65  CONTINUE 

4.  Restart  Field  For  Continuation  Runs.  (Ignored  if  not  a 
restart  case.)  Flowfield  information  for  restarting  a  computation 
is  obtained  by  reading  the  conservative  variables  data  on  TAPE7. 
This  tape  contains  the  local  time  step,  DELTAU(I,J)  and  the 
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conservative  variables,  RHO(I,J),  RHOU(I,J),  RH0V(I,J),  and  E(I,J) 
for  each  grid  point.  This  file  is  written  upon  completion  of  any 
run  and  is  then  read  from  TAPE7  by  subroutine  INITIA  for  use  in 
restarting.  The  format  used  is: 

READ ( 7 )  ( (DELTAU(I,J) ,  1=1, IL) ,  J=1,JL) 

READ (7)  ( (RHO(I, J) ,  RHOU(I,J),  RHOV(I,J), 

E (I , J) ,  1=1, IL) ,  J=l, JL) 


Output  Data  Description 

1.  Printed  Results.  The  printed  output  is  the  same  as  for 
the  perfect  gas  Euler  code  described  in  Output  Data  Description. 

2.  Namelist  INPUT.  (If  not  a  restart  case,  namelist  INPUT 
is  written  out  to  TAPE 2 . )  In  this  output,  the  value  of  IRVN  is 
changed  by  the  code  from  0  to  1  to  prepare  for  the  restart.  This 
change  enables  the  next  run  to  read  the  convergence  rate  from  TAPE9 
and  the  conservative  variables  from  TAPE7  as  determined  by  the 
previous  run.  The  conservative  variables  are  then  used  as  the 
initial  condition  for  the  succeeding  run. 

3.  Convergence  Rate.  This  output  presents  the  convergence 
rate  of  the  TLNS  equations,  SS(K),  K=l,4  for  each  iteration.  The 
convergence  rates  from  the  previous  run  are  stored  on  TAPE10  by 
subroutine  INITIA,  so  that  a  continuous  concatenation  of  the 
convergence  will  be  available.  The  format  is: 

WRITE  (10,  502)  NDUM,  (SS(K),  K=l,4) 

502  FORMAT  (15,  3X,  4(1X,  E14.7)) 

The  convergence  levels  from  the  current  run  are  added  to  this 
same  file  by  subroutine  SUPPLY  (ENTRY  CHECK) . 

4.  Conservative  Variables.  The  conservative  variables  used 
for  restart  are  written  on  TAPE 8  in  the  same  manner  described  in 
Output  Data  Description  for  the  perfect  gas  Euler  code. 

Choosing  the  Laminar  Viscosity 

The  laminar  viscosity  is  chosen  from  among  the  three  choices 
in  Eqn.  54  by  subroutine  MU  (ENTRY  MULAN  (II)).  One  of  the 
following  three  lines  of  FORTRAN  must  be  included  by  recompilation. 
For  the  constant  viscosity  case,  set, 

ZMU(J)  =  ZMUO , 

For  the  Sutherland  law  case,  set, 

ZMU(J)  =  ZMUO*TOS/TTS* (TT/TREF) **1.5, 
and  for  the  power  law  case,  set, 

ZMU(J)  =  ZMUO* (TT/TREF) **OMEGA 

Viscous  Flow  of  Real  Gases,  Code  DDADIPBG 

The  code  DDADIPBG  solves  the  TLNS  equations  for  a  real  gas 
flow.  The  molecular  weight,  Mw,  and  the  internal  energy,  e,  are 
functions  of  the  density,  and  the  temperature  as  noted  before,  and 
representative  properties  for  air  have  been  used  for  verification 
purposes  as  stated  earlier.  The  laminar  viscosity,  p,  and  the 
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laminar  thermal  conductivity,  k,  are  given  as  simplified  curve  fit 
functions  of  the  density,  p,  and  the  internal  energy,  e.  The 
turbulent  viscosity,  y+,  is  determined  by  the  Baldwin-Lomax  model 
of  turbulence9.  The  turbulent  thermal  conductivity,  k-p,  is  again 
determined  from  the  turbulent  viscosity  using  the  turbulent  Prandtl 
number. 

In  this  version  of  the  code,  the  number  of  mesh  points  allowed 
under  the  NOS-BE  system  is  limited  to  60  points  in  the  £  direction 
by  40  points  in  the  T]  direction.  The  code  is  designed  with  restart 
capability  so  that  converged  solutions  can  be  obtained  in 
sequential  runs. 

Input  Data  Description 

1.  Namelist  INPUT.  The  variables  in  this  namelist  are 
identical  to  those  for  the  perfect  gas  TLNS  code  (version  DDADIPBC, 
Input  Data  Description  for  the  perfect  gas  Euler  code)  except  that 
the  real  gas  variables,  RG,  AMWO,  and  GAMMAO  which  were  used  in  the 
real  gas  Euler  code  (version  AXI2DSG,  Input  Data  Description  for 
the  real  gas  Euler  code)  are  also  included  here. 

2.  Namelist  DINPL.  This  namelist  is  used  to  read  in  the 
real  gas  properties.  It  is  identical  to  the  namelist  used  in  the 
real  gas  Euler  code,  version  AXI2DSG,  in  Input  Data  Description  for 
the  real  gas  Euler  code. 

3.  Mesh  Point  Coordinates.  The  transformed  coordinates  are 
again  read  in  in  the  same  fashion  as  before.  The  input  is 
described  in  Input  Data  Description  for  the  perfect  gas  version  of 
the  Euler  code. 

4.  Convergence  Rate  Summary  for  Restart  Cases.  The  table 
of  convergence  levels  from  previous  runs  is  again  used  as  input  for 
restart  cases  so  a  single  concaterated  file  of  the  error  at  each 
iteration  is  available.  Input  details  are  given  in  Input  Data 
Description  for  the  perfect  gas  version  of  the  TLNS  code. 

5.  Conservative  Variables  Input  for  Restart.  Restarting 
the  code  is  accomplished  by  reading  in  the  conservative  variables 
from  the  last  time  step  of  the  previous  run.  Details  of  this 
output  file  are  given  in  Output  Data  Description  of  the  real  gas 
version  of  the  Euler  code.  Comments  as  to  re-entering  this  data 
for  restart  purposes  are  given  in  Input  Data  Description  of  the 
perfect  gas  version  of  the  TLNS  code. 

Output  Data  Description 

Output  from  the  real  gas  version  of  the  TLNS  code  is  composed 
of  four  parts.  Flowfield  results  are  available  in  either  long  or 
short  form  and  are  written  to  TAPE6  as  described  in  Input  Data 
Description  (real  gas  version  of  the  perfect  gas  code) .  A  printout 
of  Namelist  INPUT  is  written  to  TAPE9  as  described  in  Output  Data 
Description  of  the  perfect  gas  TLNS  code,  with  the  value  of  IRVN 
reset  as  necessary.  The  convergence  rate  table  including  a  summary 
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of  all  runs  to  date  is  written  to  TAPE7  for  plotting  or  restarting, 
and  the  conservative  variables  for  the  last  time  step  computed  are 
written  to  TAPE8  for  similar1  purposes.  Format  details  and  contents 
of  these  tapes  are  described  in  Output  Data  Description  of  the 
perfect  gas  Euler  code. 
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EXAMPLE  CALCULATIONS 


A  broad  spectrum  of  perfect  gas  calculations  are  given  in 
Appendix  D  along  with  some  careful  validations  of  the  code 
predictions  against  independent  standards.  The  perfect  gas 
examples  include  inviscid  calculations  for  the  transonic  flow  in 
the  converging  and  throat  sections  of  various  nozzle  geometries,  as 
well  as  inviscid  solutions  for  the  divergence,  supersonic  portion 
of  various  nozzles.  In  addition,  similar  perfect  gas  results  are 
presented  to  illustrate  TLNS  solutions  in  both  the 
subsonic-transonic  and  the  supersonic  sections  of  the  nozzle. 

These  results  include  both  laminar  and  turbulent  conditions.  In 
addition,  Appenuix  D  gives  both  inviscid  and  viscous  solutions  for 
swirling  flow  in  the  divergent  section  of  various  nozzles.  Results 
concerning  the  effect  of  back  pressure  on  the  flowfield  structure 
near  the  nozzle  lip  are  presented  for  moderate  variations  in  back 
pressure  showing,  in  particular,  the  manner  in  which  recirculating 
regions  are  set  up  near  the  nozzle  wall.  Finally,  examples  showing 
the  coupled  effect  of  nozzle  wall  cooling  on  wall  temperature  and 
heat  flux  are  presented. 

In  the  present  section,  we  summarize  some  real  gas 
calculations  for  one  particular  nozzle  geometry  and  we  present  some 
complete  subsonic-transonic-supersonic  solutions  for  an  entire 
converging-diverging  nozzle. 

Real  Gas  Calculations  For  The  R-S  Nozzle 

The  coordinates  of  the  supersonic  section  of  3  specific 
contoured  nozzle  are  given  in  Table  5.  This  nozzle  geometry  was 
used  for  the  real  gas  calculations  given  in  this  section.  The 
Namelist  INPUT  data  needed  to  run  the  cases  described  herein  is 
given  in  Figs.  1  through  4  for  the  perfect  gas  Euler  version  of  the 
code  (Fig.  1) ,  the  real  gas  Euler  version  (Fig.  2) ,  the  perfect  gas 
TLNS  version  (Fig.  3) ,  and  the  real  gas  TLNS  version  (Fig.  4) . 

These  calculations  were  run  on  the  AFAL  CYBER  180/840  machine  under 
the  NOS-BE  operating  system.  This  available  core  size  with  this 
system  is  very  limited  so  that  maximum  allowable  grid  sizes  are 
quite  sparse.  Nevertheless,  the  data  give  an  indication  of  the 
code's  validity.  The  real  gas  calculations  were  all  done  for 
equilibrium  air  as  noted  earlier. 

Figure  5  shows  the  nozzle  geometry  and  the  70  x  44  grid  that 
was  used  for  the  Euler  equation  solutions.  Figure  6  shows  Mach 
number  and  u-velocity  component  contour  plots  for  the  perfect  gas 
constant  specific  heat  calculation.  Comparisons  of  calculations 
made  with  the  perfect  gas  code  of  Appendix  D  and  the  real  gas  Euler 
code  with  property  tables  that  correspond  to  constant  specific 
heats  were  identical  indicating  the  real  gas  routines  were  working 
properly.  These  latter  perfect  gas  calculations  are  shown  in  Fig. 
7.  The  effects  of  real  gas  properties  (for  equilibrium  air)  are 
shown  on  Fig.  8.  As  would  be  expected,  real  gas  properties 
drastically  diminish  the  Mach  numbers,  but  the  exit  velocity 
profiles  are  slightly  higher  for  the  real  gas  case  partially 
because  of  the  low  value  of  Y  used  for  the  perfect  gas 
calculations. 
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Thin-layer  Navier-Stokes  calculations  are  shown  on  Figs.  9 
through  11.  Figure  9  shows  the  60  x  40  grid  used  for  these 
calculations  (refer  to  Table  2) ,  while  Fig.  10  shows  TLNS  results 
based  on  perfect  gas  properties  and  Fig.  11  shows  results  for  the 
real  gas  calculation.  Again,  the  effect  of  real  gas  properties  is 
to  reduce  the  expansion  Mach  number. 

Additional  Perfect  Gas  Calculations  For  The  Converging-Diverging 
Nozzle 

As  an  additional  calculation  demonstrating  the  coupled 
transonic-supersonic  codes,  we  have  computed  inviscid  and  viscous 
solutions  for  a  complete  converging-diverging  nozzle.  Figure  12 
shows  the  converging  portion  of  the  nozzle  and  the  120  x  80  grid 
used  for  inviscid  calculations.  Figure  13  shows  the  Mach  number 
contours  in  the  subsonic  portion  of  the  nozzle.  This  transonic 
flow  solution  was  used  to  obtain  starting  information  for  the 
supersonic  portion  of  the  grid  which  is  shown  in  Fig.  14.  In  the 
diverging  section,  we  used  a  145  x  80  grid.  Figure  15  shows  the 
Mach  number  contours  in  the  supersonic  portion  of  the  nozzle,  while 
Fig.  16  shows  the  Mach  number  contours  in  the  entire  nozzle.  The 
presence  of  a  fairly  strong  shock  in  this  nozzle  gives  rise  to 
wiggles  in  the  Mach  number  contours.  These  could  be  minimized  by 
the  addition  of  a  second  order  viscosity  in  the  vicinity  of  the 
shock  or  a  TVD  limiter.  These  have  not  been  added  because  nozzle 
flows  typically  contain  only  weak  shocks.  Corresponding  pressure 
contours  for  the  inviscid  case  are  given  on  Fig.  17.  These 
calculations  are  performed  for  a  specific  heat  ratio  Y  =  1.17. 

The  subsonic  portion  of  the  grid  used  for  the  viscous 
calculations  was  increased  to  120  x  100  to  improve  resolution  in 
the  boundary  layer.  This  viscous  grid  is  shown  on  Fig.  18,  and  the 
corresponding  flow  Mach  number  contours  are  given  on  Fig.  19. 

(Note,  these  viscous  calculations  are  all  for  the  Baldwin-Lomax 
turbulence  model . ) 

The  viscous  grid  vsed  in  the  supersonic  portion  of  the  nozzle 
included  a  145  x  100  grid  as  shown  in  Fig.  20.  The  corresponding 
Mach  number  contours  in  the  supersonic  section  are  given  in  Fig. 

21,  while  Mach  number  contours  for  the  complete  nozzle  are  given  in 
Fig.  22.  Note,  that  comparison  of  Figs.  16  and  22  shows  that  the 
inviscid  flow  expands  slightly  more  than  the  viscous  case  as  would 
be  expected.  These  calculations  are  for  a  throat  Reynolds  number 
of  5  x  105.  Corresponding  static  pressure  contours  for  the  entire 
nozzle  are  given  in  Fig.  23. 

Finally,  the  convergence  of  the  viscous  TLNS  calculation  is 
presented  in  Fig.  24.  This  figure  shows  that  the  RMS  errors  in  the 
solution  were  reduced  by  more  than  five  orders  of  magnitude  in  40 
iterations . 
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Table  1.  SUMMARY  OF  THE  VARIOUS  VERSIONS  OF  THE  CODE 


Code 

Name 

Governing 

Equations 

Gas 

Molecular 

Weight 

Mw 

Internal 

Energy 

e 

Viscosity 

Thermal 

Conduct¬ 

ivity 

k* 

AXI2DSA 

Euler 

Per¬ 

fect 

Const . 

CVT 

- 

- 

AXI2DSG 

Euler 

Real 

Mw(p,T) 

e(p,T) 

- 

- 

DDADIPBC 

TLNS 

Per¬ 

fect 

Const. 

CVT 

Power  Law 
or  Suther¬ 
land's  law 

— 

DDADIPBG 

TLNS 

Real 

MW(P,T) 

e(p,T) 

H(p,e) 

k  ( p , e) 

*Related  to  viscosity  by  Prandtl  number. 
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Table  2 .  CHARACTERISTICS  OF  CODE  EXECUTION 


Code 

Name 

Number  of 

Mesh  Points* 

Number  of 
Iterations 

Executing  Time 
on  CYBER  180/840 
Per  Iteration  Per 
Mesh  Point 
(Milliseconds) 

AXI2DSA 

70  X  44 

30 

3.851 

AXI2DSG 

70  X  44 

30 

7.412 

DDADIPBC 

70  X  44 

400 

27.211 

DDADIPBG 

60  X  44 

2000 

121.912 

*Maximum  storage  allowed  by  NOS-BE  system.  NOS-VE  allows  much 
denser  grids. 
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Table  3 .  UNITS  FOR  INPUTS 


Physical  Quantity 

Unit 

Symbol 

Length 

meter 

m 

Mass 

kilogram 

kg 

Density 

kilogram  per 
cubic  meter 

kg/m3 

Time 

second 

sec 

Temperature 

kelvin 

K 

Force 

newton 

N 

Pressure 

newton  per 
square  meter 

n/m3 

Energy 

joule 

J 

Molecular  weight 

kg  per 

kilogram  mole 

kg/kg  mole 

Viscosity 

newton-second 
per  square  meter 

(n*  sec)/m2 

Thermal 

conductivity 

joule 

kelvin-meter-second 

J/ (k*  m. sec) 

Table  4.  FUNCTION  SUBPROGRAMS  FOR  REAL  GAS  PROPERTY  EVALUATION 


Function 

Subprogram 

Name  and 
Arguments 

Property 

Computed 

i - 

Comments 

FAMW  (RHO,T) 

Mw ( P / T ) 

Molecular  weight 

FE  (RHO,  T) 

E ( p , T) 

Internal  energy 

FT  (RHO,  E) 

T(p,e) 

Temperature 

FDMDRT  (RHO,  T) 

(3mw/3P)t 

Molecular  weight  derivative  as 
function  of  p  and  T 

FDMDTR  (RHO,  T) 

(3mw/3t) 

P 

Molecular  weight  derivative  as 
function  of  p  and  T 

FDTDRE  (RHO,  E) 

(3T/9p)e 

Temperature  derivative;  function 
of  p  and  e 

FDTDER  (RHO,  E) 

(3T/3e) 

P 

Temperature  derivative;  function 
of  p  and  e 

FZMU  (RHO,  E) 

V(p,e) 

Viscosity 

FZK  (RHO,  E) 

k(p,e) 

Thermal  conductivity;  computed 
from  Prandtl  number 

FDMUDRE  (RHO,  E) 

(3p/3p)e 

Viscosity  derivative;  function  of 
p  and  e 

FDMUDER  (RHO,  E) 

(3p/3e) 

P 

Viscosity  derivative;  function  of 
p  and  e 

FDKDRE  (RHO,  E) 

(3k/3p) e 

Thermal  conductivity  derivative 

FDKDER  (RHO,  E) 

(3k/3e) 

P 

Thermal  conductivity  derivative 

Table  5a.  OCWIOJR  POINTS  OF  THE  RS  NOZZLE  (70.  441 


I 

YC 

XC 

ALFHA 

I 

YC 

XC 

Ai_IHA 

1 

.01006 

0.00000 

0.00000 

41 

. 07105 

. 11669 

21.64860 

2 

.01009 

.00093 

3.36523 

42 

.07377 

.12362 

21.08610 

3 

.01018 

.00191 

7.10866 

43 

.07654 

.13090 

20.52531 

4 

.01035 

.00295 

11.07978 

44 

.07936 

. 13857 

19.96692 

5 

.01061 

.00403 

15.31548 

45 

. 08224 

.14662 

19.41090 

6 

.01097 

.00517 

19.86431 

46 

.08518 

.15509 

18.85670 

7 

.01147 

.00637 

24.79635 

47 

.08818 

.16400 

18.30537 

8 

.01213 

.00763 

30.20430 

48 

.09123 

.17336 

17.75549 

9 

.01301 

.00896 

36.34891 

49 

.09433 

.18321 

17.20873 

10 

.01403 

.01035 

36.33356 

50 

.09748 

.19356 

16.66520 

11 

.01511 

.01182 

36.27180 

51 

.10068 

.20444 

16.12458 

12 

.01623 

.01336 

36.11491 

52 

.10393 

.21588 

15.58597 

13 

.01741 

.01498 

35.89578 

53 

.10723 

.22792 

15.05187 

14 

. 01864 

.01668 

35.62613 

54 

.11056 

.24056 

14.52032 

15 

.01992 

.01847 

35.30239 

55 

.11395 

.25386 

13.99123 

16 

.02124 

.02035 

34.93930 

56 

.11736 

.26785 

13.46773 

17 

.02261 

.02233 

34.54235 

57 

.12081 

.28255 

12.94699 

18 

.02404 

.02442 

34.11646 

58 

.12429 

.29801 

12.42823 

19 

.02551 

.02660 

33.66482 

59 

.12780 

.31426 

11.91538 

20 

.02703 

.02890 

33.19136 

60 

.13133 

.33135 

11.40440 

21 

.02859 

.03132 

32.70377 

61 

.13487 

.34931 

10.89921 

22 

.03021 

.03387 

32.20513 

62 

.13841 

.36820 

10.39807 

23 

.03188 

.03654 

31.69175 

63 

.14197 

.38806 

9.89894 

24 

.03360 

.03935 

31.16775 

64 

.14552 

.40894 

9.40745 

25 

.03537 

.04231 

30.63482 

65 

.14907 

.43090 

8.91692 

26 

.03719 

.04542 

30.09359 

66 

.15258 

.45398 

8.43359 

27 

.03906 

.04869 

29.54643 

67 

.15609 

.47825 

7.95193 

28 

.04099 

.05212 

28.99329 

68 

.15957 

.50377 

7.48025 

29 

.04297 

.05573 

28.43651 

69 

.16286 

.53060 

7.00226 

30 

.04500 

.05953 

27.87618 

70 

.16596 

.55880 

6.52500 

31 

.04709 

.06353 

27.31309 

32 

.04923 

.06772 

26.74888 

33 

.05143 

.07214 

26.18210 

34 

.05368 

.07678 

25.61507 

35 

.05599 

.08166 

25.04"16 

36 

.05836 

.08679 

24.47881 

37 

. 06078 

.09219 

23.91089 

38 

.06325 

.09786 

23.34445 

39 

.06580 

.10382 

22.77762 

40 

.06840 

.11009 

22.21202 

Note:  XC  and  YC  are  axial  and  radial  coordinates  in  meters; 
ALFHA  is  wall  angle  in  degrees. 
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Table  5b.  CONTOUR  POINTS  OF  TOE  RS  NOZZIE  (60.  4CH 


T 

YC 

XC 

ALPHA 

I 

YC 

XC 

ALPHA 

1 

.01006 

0.00000 

0.00000 

41 

.08717 

.16098 

18.48817 

2 

.01009 

.00093 

3.36523 

42 

.09087 

.17226 

17.81855 

3 

.01018 

.00193 

7.15475 

43 

.09466 

.18427 

17.15137 

4 

.01036 

.00298 

11.22513 

44 

.09852 

.19705 

16.48805 

5 

.01063 

.00411 

15.62255 

45 

.10246 

.21065 

15.82853 

6 

.01102 

.00531 

20.41094 

46 

.10647 

.22513 

15.17335 

7 

.01157 

.00658 

25.67889 

47 

.11056 

.24053 

14.52154 

8 

.01232 

. 00794 

31.56032 

48 

.11470 

.25693 

13.87519 

9 

.01332 

.00938 

36.34023 

49 

.11891 

.27438 

13.23191 

10 

.01445 

.01092 

36.32034 

50 

.12317 

.29295 

12.59448 

11 

.01565 

.01255 

36.20025 

51 

.12747 

.31272 

11.96280 

12 

.01692 

.01429 

35.99863 

52 

.13181 

.33376 

11.33653 

13 

.01826 

.01615 

35.71517 

53 

.13617 

. 35615 

10.71358 

14 

.01966 

.01812 

35.36682 

54 

.14055 

. 37998 

10.09854 

15 

.02114 

.02021 

34.96643 

55 

.14492 

.40534 

9.49067 

16 

.02269 

.02245 

34.51934 

56 

.14929 

.43233 

8.88584 

17 

.02431 

.02482 

34.03213 

57 

.15363 

.46106 

8.29007 

18 

.02600 

.02735 

33.51066 

58 

.15796 

.49163 

7.70414 

19 

.02777 

.03004 

32.95256 

59 

.16210 

.52417 

7.11456 

20 

.02961 

.03291 

32.39225 

60 

.16596 

.55880 

6.52500 

21 

.03152 

.03596 

31.80246 

22 

.03351 

.03920 

31.19516 

23 

.03557 

.04266 

30.57319 

24 

.03772 

.04633 

29.93819 

25 

.03994 

.05025 

29.29224 

26 

.04225 

.05441 

28.63803 

27 

.04463 

.05884 

27.97653 

28 

.04710 

.06356 

27.30897 

29 

.04966 

.06858 

26.63766 

30 

.05230 

.07392 

25.96185 

31 

.05502 

.07960 

25.28350 

32 

.05784 

.08565 

24.60259 

33 

.06074 

.09209 

23.92056 

34 

.06373 

.09895 

23.23896 

35 

.06681 

.10624 

22.55594 

36 

.06998 

.11400 

21.87451 

37 

.07324 

.12227 

21.19345 

38 

.07659 

.13106 

20.51365 

39 

.08003 

.14042 

19.83620 

40 

.06356 

.15038 

19.16112 

Note:  XC  and  YC  are  axial  and  radial  coordinates  in  meters; 
ALPHA  is  wall  angle  in  degrees. 
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I 


-$INFUr 


IL 

=  70, 

JL 

=  44, 

ITRl 

=  1, 

ITRN 

=  1, 

ITRS 

=  1, 

ISUP 

=  o, 

CFL 

=  5., 

ITIME 

=  1, 

THETA 

=  1., 

RMl 

=  1.02, 

RM2 

=  8., 

AIN 

=  0.05, 

AEX 

=  0.275, 

RL 

=  0.8485, 

PO 

=  3509431.4 

TO 

=  3497.8, 

CP 

=  1308.79, 

CV 

=  1021.77, 

GAMMA 

=  1.2809, 

GM1 

=  .2809, 

INER 

=  30, 

IB1 

=  2, 

IB2 

«  70, 

NORD 

=  1, 

CMEGAX 

=  0.0, 

CMEGAY 

=  .5, 

IREAD 

=  1, 

IWRT 

=  1, 

IRJN 

=  1, 

$END 

Fig.  1  Sample  Input  For  AXI2DSA:  Namelist  INPUT. 
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-$INR7T 


IL 

=  70, 

JL 

=  44, 

ITR1 

=  1, 

FERN 

=  1, 

ITRS 

=  1, 

ISUP 

=  o, 

CFL 

=  5., 

itime 

=  1, 

THETA 

-  1., 

EMI 

=  1.02, 

EM2 

=  8., 

AIN 

=  0.01, 

AEX 

=  0.11, 

RL 

*  0.245, 

PO 

=  3509431.4 

TO 

=  3497.8, 

RG 

=  8314., 

AMMO 

=  28.9670, 

GAMMftO 

=  1.2809, 

INER 

=  30, 

IB1 

*  2, 

IB2 

=  70, 

NORD 

=  1, 

CMEGAX 

=  0.0, 

CMBGAY 

=  .5, 

IREAD 

“  1/ 

IWRT 

-  1, 

IRUN 

$END 

-  1, 

Namelist  INHJT 


-$INRJT 


IL 

= 

70, 

JL 

= 

44, 

NBBG 

= 

1, 

NEND 

1, 

NITER 

= 

100, 

THETA 

= 

1.0, 

NORD 

1, 

CFL 

= 

4.0, 

CELL 

= 

1.0, 

ITIME 

— 

1, 

CMEGAX 

= 

0.0, 

CMEGAY 

= 

0.5, 

AIN 

= 

0.05, 

AEX 

= 

0.275, 

RL 

= 

0.8485, 

EST 

= 

0.0, 

FSTY 

= 

0.0, 

EMI 

= 

1.02, 

RM2 

= 

8.0, 

IVISC 

= 

1, 

IWALL 

= 

o, 

GAMMA 

= 

1.2809, 

CP 

= 

1308.79, 

REN 

= 

1.0, 

FRN 

= 

0.7, 

FRET 

= 

0.9, 

TREF 

= 

3109.55, 

73/VO 

= 

1.02319E-05 

CMEGA 

= 

0.8130, 

PO 

= 

3509431.4, 

TO 

= 

3497.8, 

TWAIjL 

= 

0.0, 

PB 

= 

0.0, 

IREAD 

= 

1, 

IWRT 

0, 

IEUN 

$END 

0, 

jg-jiL  «J  ll1 


>r  DDADIFBC:  Namelist  TNH7T 


-$iNKrr 


IL 

=  60, 

JL 

=  40, 

NBBG 

=  1, 

NEND 

=  100, 

NITER 

=  100, 

THETA 

=  1.0, 

N0RD 

=  0, 

CFL 

=  0.5, 

CFLl 

=  1.0, 

ITIME 

=  0, 

CMEGAX 

=  0.0, 

CMBGAY 

=  0.08, 

AIN 

=  . 5E-01 , 

AEX 

=  . 275E+00, 

RL 

=  . 8485E+00 , 

FST 

=  0.0, 

FSTY 

=  0.0, 

FM1 

=  1.02, 

RM2 

=  5.0, 

IVISC 

=  1, 

IWALL 

=  o, 

RG 

=  8314., 

AMWO 

=  28.967, 

GAMMAO 

=  1.2809, 

CP 

=  1308.79, 

REN 

=  . 1E+01 , 

ERN 

=  0.7, 

FROT 

=  0.9, 

TREF 

=  1.0, 

ZMJO 

=-  .  102319E-I 

CMEGA 

=  .838551, 

PO 

=  3509431.4 

TO 

=  3497.8, 

TVJALL 

=  0.0, 

PB 

=  0.0, 

TREAD 

=  1, 

IWRT 

=  o. 

IHJN 

=  o, 

$END 

Fia.  4  Sairole  Input  For  DDADIFBG;  Namelist  INPUT. 
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Mach  Number  Contours 


u -  Velocity  Contours 


Fig.  6  Contours  of  Mach  Number  (Top)  and  U-Velgcity_( 
Perfect  Gas  Solution  of  Euler  Equations. 


u -Velocity  Contours 


Fig. 


Contours  of  Mach  Number  (Top)  and  U-Velocitv  (Bottom) 
Real  Gas  Calculations  of  Euler  Equations  Based  on 
Equilibrium  Air. 
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For 


Fig.  9  Nozzle  Geometry  and  60  x  40  Grid  Used  for  Thin-Laver 


Navier-Stokes  Calculations. 


1.2  2.0  3.0 


Kav.li  Number  Contours 


3-0 


*4.0 


1.0  2.0 


Fia.  22  Mach  Number  Contours  of  Viscous  Calculations. 
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APPENDIX  A 


JACOBIAN  MATRICES  AND  ASSOCIATED  ALGEBRA  FOR  REAL  GAS  EQUATIONS 

The  present  Appendix  summarizes  a  number  of  important 
algebraic  quantities  for  the  real  gas  equations.  In  all  cases,  we 
indicate  how  these  quantities  simplify  for  the  perfect  gas  case. 

We  start  by  specifying  the  Jacobian  matrices  of  the  inviscid  flux 
vectors  E  and  F.  These  are  defined  right  below  Eqn.  1  in  the  main 
body  of  the  report.  These  Jacobians  are  respectively,  A  =  9e/9q 
and  B  =  9f/9q.  They  are  composed  of  the  following  terms: 


K 


K 


A  or  B  = 


3p 

-UU+KX  f 


-VU+K 


9p 

3p 


u+Vu+3T?^j- 

Kxv+Ky  aw 


Kyu+Kx  tj(pv) 

U+Vv+  stW1 


K 


9p 

3T 


K  9p_ 

Ky  9e 


U( 


3p  en+P 
3p"  P 


) 


u 


3  (pu)  Kx  p 


T!  dP  ik 

U9  (pv)  Ky  p 


O 

(A.l) 


where  U  =  KXU  +  KyV,  and  K  =  £,  U  =  U  for  A,  while  K  =  r\ ,  U  =  V  for 
B. 

For  the  real  gas  case,  the  derivatives  of  pressure  are  given 

by: 


3P  =  A  -  A 
<*P  AP  e  p(p 


2  2, 

-  u  -  v  ) 


9p 
3  pv 

9p 

3p\ 

9p 


- -  =  _  A  H 

9pu  e  p 

dpv  e  p 


=  Ae  P 


(A.2) 


where, 


P  P 

Ap  =  p  BP  +  T  BT 


A  = 
e 


”  7  bt  ^  p 


,3t 

(3p}e 


(A. 3) 


and, 


bp  ■ 


_  £_  /3mw. 

M  'cJp'T 
w  r 


bt  -  1 


mT 

W 


1 


(A. 4) 


For  the  perfect  gas  case,  Mw  is  a  constant,  so  Bp  =  Bj  =  1. 
In  addition,  (3T/3p)e  =  0,  while  (3T/3e)p  =  1/CV.  Thus,  for  the 
perfect  gas  case,  we  have, 


and, 


Ap  =  p/p 

Ae  =  (Y-l)p 


(A. 5) 


(A. 6) 


Substituting  these  values  into  the  expressions  for  the  derivatives 
of  pressure  and  inserting  in  the  matrices  A  and  B  gives  the 
standard,  perfect  gas  form  of  these  Jacobians. 

The  modal  matrices  T-1P~1  and  Tb-1p-1  are:  t^P-1  = 
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The  inverses  of  these  matrices,  PT^  and  PT3  are  given  by,  PT  = 
1 


u 


0 

Ab 

p 

rt c 

Kjp 

JAi  +  Kx> 

A<1- 

-KxP 

A<i ♦  Ky> 

1 

>1 0 
* 

Kj) 


a 


P  (KyU+K^v)  /To  ai£f.||o^(Klu+Kjv) 

> 

(A. 8) 


1 

5p 


where, 
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Again,  the  derivatives  of  pressure  are  as  given  above,  and  the 
limiting  values  for  perfect  gases  again  give  the  familiar  results 
of  Appendix  C. 

For  the  right-hand  side  terms,  we  need  the  Jacobians  D'  and 

D' '  . 

For  the  source  term,  the  Jacobian  of  the  algebraic 
(non-differentiated)  terms  in  H  is, 
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where,  again  the  pressure  derivatives  are  given  above  for  both  real 
and  perfect  gas  variables. 

Similarly,  the  Jacobian  of  the  differentiated  (first  order) 
terms  in  H  is, 
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where , 
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The  derivatives  of  y  (and  later  k  and  T)  with  respect  to  Q 
require  some  comment.  By  taking  y  =  y(p,e)  and  k  =  k(p,e)  and 
using  the  chain  rule,  we  obtain, 

-  <g“)e  dp  +  <g“>p  de  (A. 13) 


with  a  similar  expression  for  the  changes  in  k  and  T.  Derivatives 
of  y ,  k,  or  T  are  then  defined  in  terms  of  the  derivatives  (3y/3p)e 
and  (3y/3e)p  etc.  which  are  obtained  from  the  property  tables. 
(Note,  the  dy/3p)  appearing,  for  example,  in  Eqn.  A. 12a  is  not 
equal  to  (3y/3p)e  because  the  term  in  Eqn.  A. 12a  is  obtained  with  Q 
as  the  independent  variable,  whereas  the  latter  considers  p  and  e 
as  independent  variables. 
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while  the  matrix  B^  is  given  by. 
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The  derivatives  of  temperature  are  first  expanded  by  the  chain  rule 
expression  given  in  Eqn.  A. 13  above.  Because  temperature  is  a 
primary  variable  of  interest,  we  re-express  this  here  as, 
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and  we  introduce  the  coefficients  Cp  and  ce  such  that, 
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In  terms  of  these  variables,  we  can  express  the  four  derivatives  of 
temperature  as: 


9t 

3p  = 


c 


p 


(-  |°  +  (u2  +  v2)) 

(A. 17a) 

0T 

d(pu) 

u 

"  Ce  p 

(A. 17b) 

9t 

d(pv) 

(A. 17c) 

9t 

1 

(A. 17d) 

JT  = 

n 

Ce  p 

For  completeness,  we  also  define  the  speed  of  sound.  The 
speed  of  sound,  c,  is  given  by: 
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Roe  averaging  is  used  to  ensure  conservation  even  when  some 
terms  are  expressed  in  non-conservative  form.  Details  of  the 
procedure  are  given  in  Ref.  4.  In  brief.  Roe  averaging  is  given 
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The  total  internal  energy,  eQ,  is  related  to  the  stagnation 
enthalpy  by. 


phQ  =  eQ  +  p 
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where  p  is  the  pressure.  P  can  be  evaluated  from, 

P  =  1-^.1  (pho  -  i/2p(U2  +  v2)  (B.  6) 

and  assuming  an  arithmetic  averaging  for  which  is  given  by, 
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